!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!> \brief
!! Raising bubble test case.
!!
!! \n
!!
!! This test case uses a neutrally stratified atmosphere, i.e. an
!! atmosphere for which \f$\theta=const\f$ and \f$N=0\f$. The
!! hydrostatic balance yields:
!! \f{equation}{
!!  \theta(z)=\theta_s, \qquad \pi(z) = 1-\frac{g}{c_p\theta_s}z,
!! \f}
!! from which any other thermodynamical variable can be recovered.
!!
!! \bug Add details and references.
!<----------------------------------------------------------------------
module mod_warmbubble1_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_physical_constants, only: &
   mod_physical_constants_initialized, &
   t_phc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_warmbubble1_test_constructor, &
   mod_warmbubble1_test_destructor,  &
   mod_warmbubble1_test_initialized, &
   test_name,        &
   test_description, &
   phc,              &
   ntrcs,            &
   ref_prof, t_ref,  &
   coeff_init,       &
   coeff_outref

 private

!-----------------------------------------------------------------------

 ! public interface
 character(len=*), parameter   :: test_name = "warm-bubble-1"
 character(len=100), protected :: test_description(3)
 type(t_phc), save, protected  :: phc
 integer, parameter            :: ntrcs = 0
 character(len=*), parameter   :: ref_prof = "isothermal"
 real(wp), target              :: t_ref = 300.0_wp

 logical, protected :: mod_warmbubble1_test_initialized = .false.

 ! private members
 real(wp) :: &
   nu,     & ! viscosity
   dtheta, & ! potential temperature anomaly
   xc(3),  & ! bubble center
   rc_i,   & ! internal bubble radius
   sigma,  & ! bubble profile
   theta_s   ! background theta
 character(len=*), parameter :: &
   test_input_file_name = 'warmbubble1_test.in'
 character(len=*), parameter :: &
   this_mod_name = 'mod_warmbubble1_test'

 ! Input namelist
 namelist /input/ &
   nu, dtheta, xc, rc_i, sigma, theta_s

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_warmbubble1_test_constructor(d)
  integer, intent(in) :: d

  integer :: fu, ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)   .or. &
       (mod_kinds_initialized.eqv..false.)      .or. &
       (mod_fu_manager_initialized.eqv..false.) .or. &
       (mod_physical_constants_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_warmbubble1_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1)') &
     'Warm test case with smooth bubble, dimension d = ',d
   write(test_description(2),'(a,e9.3,a)') &
     '  viscosity nu = ',nu,';'
   write(test_description(3),'(a,e9.3,a)') &
     '  potential temperature anomaly d\theta = ',dtheta,'.'

   ! All the physical constants use the default values: phc unchanged

   t_ref = 300.0_wp

   mod_warmbubble1_test_initialized = .true.
 end subroutine mod_warmbubble1_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_warmbubble1_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_warmbubble1_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_warmbubble1_test_initialized = .false.
 end subroutine mod_warmbubble1_test_destructor

!-----------------------------------------------------------------------

 !> Potential temperature perturbation with hydrostatic balance.
 !!
 !! The reference state u_r at the same points x where the initial
 !! condition is desired must be passed. In fact, once prescribed the
 !! potential temperature perturbation, the perturbations in the
 !! conservation variables can not be computed without the reference
 !! state, due to the nonlinearity of the problem.
 !<
 pure function coeff_init(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1),size(x,2))

  integer :: l
  real(wp) :: r, dthetal
  real(wp), dimension(size(x,2)) :: z, theta, pi, p, t

   ! first define a neutral atmosphere: \theta = const
   z = x(size(x,1),:) ! vertical coordinate
   theta = theta_s
   pi = 1.0_wp - phc%gravity/(phc%cp*theta_s) * z

   ! now perturb the potential temperature (but leave pi unchanged)
   do l=1,size(x,2)
     r = sqrt(sum((x(:,l)-xc(1:size(x,1)))**2))
     if (r < rc_i) then
       dthetal = dtheta
     else
       dthetal = exp(-((r-rc_i)/sigma)**2)*dtheta
     end if
     theta(l) = theta(l) + dthetal
   enddo

   ! translate into primitive variables
   p = phc%p_s*pi**(1.0_wp/phc%kappa)
   t = pi*theta

   ! finally translate the result in conservative variables
   u0(1,:) = p/(phc%rgas*t)                     ! rho
   u0(2,:) = u0(1,:)*(phc%cv*t + phc%gravity*z) ! e
   u0(3:,:) = 0.0_wp                            ! U
   ! subtract the reference state
   u0(1,:) = u0(1,:) - u_r(1,:) ! rho
   u0(2,:) = u0(2,:) - u_r(2,:) ! e

 end function coeff_init
 
!-----------------------------------------------------------------------

 !> Neutral stability atmosphere
 pure function coeff_outref(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1),size(x,2))

  real(wp), dimension(size(x,2)) :: z, theta, pi, p, t

   ! first define a neutral atmosphere: \theta = const
   z = x(size(x,1),:) ! vertical coordinate
   theta = theta_s
   pi = 1.0_wp - phc%gravity/(phc%cp*theta_s) * z

   ! translate into primitive variables
   p = phc%p_s*pi**(1.0_wp/phc%kappa)
   t = pi*theta

   ! finally translate the result in conservative variables
   u0(1,:) = p/(phc%rgas*t)                     ! rho
   u0(2,:) = u0(1,:)*(phc%cv*t + phc%gravity*z) ! e
   u0(3:,:) = 0.0_wp                            ! U
   ! subtract the reference state
   u0(1,:) = u0(1,:) - u_r(1,:) ! rho
   u0(2,:) = u0(2,:) - u_r(2,:) ! e

 end function coeff_outref

!-----------------------------------------------------------------------

end module mod_warmbubble1_test

